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Abstract— The intent of this paper is to discuss the history and origins of Lagrangian hydrodynamic methods 
for simulating shock driven flows. The majority of the pioneering research occurred within the Manhattan 
Project. A range of Lagrangian hydrodynamic schemes were created between 1943 and 1948 by John von 
Neumann, Rudolf Peierls, Tony Skyrme, and Robert Richtmyer. These schemes varied significantly from 
each other; however, they all used a staggered-grid and finite difference approximations of the derivatives in 
the governing equations, where the first scheme was by von Neumann. These ground-breaking schemes were 
principally published in Los Alamos laboratory reports that were eventually declassified many decades after 
authorship, which motivates us to document the work and describe the accompanying history in a paper that 
is accessible to the broader scientific community. Furthermore, we seek to correct historical omissions on the 
pivotal contributions made by Peierls and Skyrme to creating robust Lagrangian hydrodynamic methods for 
simulating shock driven flows. Understanding the history of Lagrangian hydrodynamic methods can help 


explain the origins of many modern schemes and may inspire the pursuit of new schemes. 


l. PROLOGUE 


We will start by providing an overview of the history 
surrounding the hydrodynamic methods work between 
1943 and 1948 to establish context and communicate the 
sequence of events that occurred over that time frame to 
yield a suite of Lagrangian hydrodynamic methods. Sub- 
sequent sections in the paper will provide further details 
on the governing formulations and the four distinct numer- 
ical methods by John von Neumann, Rudolf Peierls, Tony 
H.R Skyrme, and then Robert Richtmyer. 

Researchers at Los Alamos seemed surprisingly ready 
in early 1944 to numerically solve the partial differential 
equations for shock hydrodynamics. Looking back to 
1940 through 1943, this readiness is not surprising. With 
the advent of World War II (WWII) the effects of the blast 
from bombs and artillery became of intense interest. Blast 
waves are described by a set of non-linear, time-dependent, 
partial differential equations. Further, the shock wave in 
an inviscid fluid has a discontinuous solution. These equa- 
tions can only be solved analytically for a few very special 
cases, in general the solution requires numerical methods. 
Exploration of blast wave effects were organized in both 
the United States (U.S.) and the United Kingdom (U.K.). 
In the U.S., the National Defense Research Committee 
(NDRC) was established in June 1940. Section B con- 
cerned “Bombs, Fuels, Gases, and Chemical Problems".! 
When the NDRC was superseded by the Office of Scien- 
tific Research and Development (OSRD) in June 1941, 
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parts of the NDRC Section B became Division 2, “Effects 
of Impact and Explosion", and Division 8, “Explosives”.! 
Vannevar Bush was the chairman of OSRD with James B. 


Conant as chairman of NDRC. 


In May 1942 OSRD Division 8 issued a report on the 
"general conditions for the existence of shock waves" by 
Hans Bethe, endorsed by George B. Kistiakowsky, Chair- 
man of Section B-1.? Hans Bethe was also exploring the 
pressure wave generated by underwater explosions.*^ von 
Neumann published a series of OSRD reports on high ex- 
plosive detonation waves and shock waves in 1942 and 
1943.55 See the Division 2 summary for the many reports 
on explosion effects.? 


Meanwhile in England, William G. Penney was car- 
rying out research on shock waves due to explosions, 
summarized in 1950,!? and also on the effects of blast 
waves during the London Blitz,!! Geoffrey I. Taylor was a 
member of the British delegation with long experience in 
hydrodynamics, shock waves, and explosives.!!-!^ Peierls 
was also working on the effects of blasts, correspond- 
ing with Taylor and Kistiakowsky.!! Peierls brought his 
assistants Klaus Fuchs and Skyrme to Los Alamos as 
members of the British delegation. Skyrme would make 
significant discoveries that enabled robust and accurate 
hydrodynamic calculations with strong shocks. By the 
time Project Y, the Los Alamos branch of the Manhattan 
Project, was formed in 1943 there was broad experience 
with shock hydrodynamics by the senior staff of both the 
U.S. and U.K. 


During the summer of 1943 Robert Serber and Edward 
Teller derived an analytic approach to solve the shock 
hydrodynamics problem.'? A set of IBM Punch Card Ac- 
counting Machines (PCAM) had been ordered in January 
1944 for use by Stanley Frankel and Eldred Nelson to 
solve the gun bomb neutronics problem. In March 1944 
the tasking of Frankel and Nelson was changed to numeri- 
cally solve a range of problems, including hydrodynamics 
problems. See the companion ANS article for a discussion 
of the PCAM and the computing facility. !° 

In March 1944, Hans Bethe reorganized T-Division, 
with the IBM calculations in Group T-2 under Serber.'” 


* T-1, “Hydrodynamics of Implosion, Super", Edward 
Teller 

* T-2, “Diffusion Theory, IBM Calculations, and Ex- 
periments”, Robert Serber 

* T-3, "Experiments, Efficiency Calculations, Radia- 
tion Hydrodynamics", Victor Weisskopf 

* T-4, “Diffusion Problems", Richard Feynman 

* T-5, “Computations”, Donald Flanders 


Peierls moved to Los Alamos in June 1944, replacing 
Teller as T-1 Group Leader. Skyrme joined Los Alamos in 
T-1 in July 1944. The stage was now set for researchers at 
Los Alamos to create hydrodynamic methods to simulate 
problems with strong shocks. 


ll. HISTORY 


After von Neumann visited Los Alamos in September 
and October 1943, he studied how to numerically solve the 
shock hydrodynamics equations. The first 1D Lagrangian 
hydrodynamic method was a staggered-grid finite differ- 
ence method, created by von Neumann, and was published 
in an OSRD report, AMP 108.1R, in March 1944. This 
first scheme was for 1D Cartesian coordinates and has 
not been well discussed in the broader literature, yet it 
is the foundation that most Lagrangian hydrodynamic 
methods are built upon. Soon thereafter von Neumann 
extended his method to 1D spherical coordinates, but that 
scheme was never published beyond Los Alamos labo- 
ratory reports. Von Neumann tested the method on a 
radially-outward propagating blast wave using the IBM 
PCAM at the Army Ballistics Research Laboratory (BRL) 
at Aberdeen Maryland in early 1944.5 Von Nuemann's 
early finite difference scheme suffered from oscillations 
at the shock front. At the time, there was much debate 
over the physical correctness and source of the oscilla- 
tions.'^?? Several alternate approaches would later be 
pursued to address the oscillations. Section III. describes 
von Neumann's previously unpublished spherical coordi- 
nate hydrodynamic method. 

Meanwhile Peierls visited Los Alamos in February 
1944.?! In his autobiography, Peierls recounts how he 
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had tried to solve the blast wave equation numerically: “I 
tried a step-by-step method, in which you take the initial 
state as known at a chain of points a certain distance apart, 
and then determine approximately the state of affairs a 
certain time interval later. One can then repeat the pro- 
cess".!! Peierls described his attempt to numerically solve 
the blast wave equations using the step-by-step method to 
the Los Alamos staff. The Los Alamos theoretical staff 
had been trying to solve the equations analytically.? On 
February 14, 1944 J. Robert Oppenheimer wrote to Gen- 
eral Groves about Peierls’ visit? He noted that Peierls 
had explained the “British methods” for integrating the 
blast wave equations in great detail. Oppenheimer was 
“not satisfied with the adequacy of the methods now in use 
in this laboratory” for solving the shock hydrodynamics 
equations. 


Peierls clearly did not know about the Courant stability 
condition,” he re-discovered it empirically.!! The Los 
Alamos staff also did not know of the Courant stability 
condition, so at first they referred to it as the "Peierls" 
stability condition. Von Neumann clearly did know of the 
Courant stability condition, he discussed it in his March 
1944 paper.** By 1945 Los Alamos was referring to it as 
the Courant condition. 


Peierls sent a letter to Oppenheimer on March 15, 
194425 that described “a new form of the shock wave 
equations which was suggested by Skyrme", see appendix 
A to read the original letter. These equations are simpler 
than the ones Peierls described during his February 1944 
visit to Los Alamos.” Peierls approach required knowl- 
edge about the material conditions ahead of the shock, 
which Skyrme's method did not. 


Peierls and Skyrme/??6 proposed coupling internal 


boundary conditions to the underlying hydrodynamic 
method to remove oscillations at a shock front. With 
these shock-fitting methods, a higher-order central dif- 
ference Lagrangian hydrodynamic method were used on 
either side of the shock, and then the discontinuous so- 
lutions on each side of the shock were coupled together 
with Rankine-Hugoniot jump conditions. The first shock- 
fiting Lagrangian hydrodynamic method was created by 
Peierls” and is briefly described by Skyrme!” and alluded 
to in a letter that he sent von Neumann in March 1944.7? 
The Peierls scheme was viewed unfavorably at the time. 
Skyrme says, “The main disadvantage of this method ... 
was that a cycle (n+1) could not be completed until the 
positions of the points (i-1) etc., behind the shock had 


been found by numerical integration". ? 


Skyrme created an alternative shock-fitting Lagrangian 
hydrodynamic method that accurately captured a discon- 
tinuous solution at a shock and evolved it across the mesh. 
The method by Skyrme required an iterative solver to find 


the shock states. That scheme appeared to have similar- 
ities to the one by Peierls, but the details on the Peierls’ 
scheme are limited in the available records. It is impor- 
tant to recognize that the shock fitting methods by Peierls 
and Skyrme were quite different from and predate the 
artificial viscosity based method by von Neumann and 
Richtmyer.*+7’ By mid-March 1944, researchers at Los 
Alamos had at least four approaches available to solve the 
shock hydrodynamics equations: schemes by Serber, von 
Neumann, Peierls, and Skyrme. 


Incorporating Rankine-Hugoniot jump conditions into 
a numerical scheme was revolutionary for that time pe- 
riod. Beyond the work by Skyrme and Peierls, there 
was just one other shock fitting scheme at that time. 
Howard Emmons?*?? published in 1944 a method for 
steady transonic flows over surfaces that arise in aeronau- 
tic applications. Emmons created an Eulerian approach 
that combined a stream-function based method for the 
subsonic region with a separate method for supersonic 
region. The method by Emmons differed significantly 
from the Lagrangian method by Skyrme and the one by 
Peierls that were designed for transient shock dynamics. 
Shock-fitting schemes similar to Peierls’ and Skryme’s 
were not published until many decades later, for example 
see the work by Maurizio Di Giacinto and Mauro Valo- 
rani,” the shock-fitting method implemented in the Los 
Alamos FRONTIER code;?!? the recent shock fitting 
book by Onofri and Paciorri,** or the modern moving 
Discontinuous Galerkin method with interface condition 
enforcement by Corrigan et. al??? The schemes by 
Peierls and Skyrme were never published in a journal so 
they remained somewhat lost to time and unknown to the 
broader scientific community. The details on Skyrme's 
shock-fitting scheme were documented in multiple labo- 
ratory reports and letters during and after the Manhattan 
Project. These were declassified many decades after au- 
thorship and have since been made available to the public. 
Details on Peierls’ shock-fitting method have not been 
found in the available records. Skyrme's shock fitting 
scheme is described in Section V.. 


Frankel, Nelson, and Naomi Livesay developed a hy- 
drodynamics computer program at Los Alamos for the 
soon-to-arrive IBM PCAM in March 1944. The program 
was tested in March 1944 on the Marchant calculators with 
the assistance of Richard Feynman and Nicholas Metropo- 
lis.??" The hydrodynamics code being discussed had 
no name, it was just referred to as the IBM code or IBM 
problems. 

The IBM PCAM arrived in Los Alamos on April 4, 
1944. A hydrodynamic calculation was started on the IBM 
PCAM “within a week of arrival" in early April 1944.7138 
See the companion ANS article on the PCAM and the 
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computing facility.' The calculation used von Neumann's 
spherical coordinate hydrodynamic method. Note that this 
method is not the method normally associated with von 
Neumann. The method that von Neumann and Richtmyer 
published in 1950 is based on artificial viscosity, and is 
quite different from the original method by von Neumann 
that was used in 1944.74 


The internal boundary conditions with Skryme's ap- 
proach was too complicated for the IBM PCAM and had 
to be performed by hand using Marchant calculators, with 
the results inserted into the IBM hydrodynamic calcula- 
tion at the appropriate position and time. This motivated 
research into alternative Lagrangian hydrodynamic meth- 
ods. 


As an alternative to a shock-fitting approach, Peierls 
suggested to von Neumann in a March 28, 1944 letter that 
he should add a dissipative mechanism to his hydrody- 
namic scheme? to remove the spurious oscillations near 
a shock. Peierls said, “Incidentally, have you thought at 
all about the following alternative way of avoiding discon- 
tinuity. In actual fact, the shock front has a finite width 
because of the viscosity and thermal conductivity of the 
medium. But artificially, assuming a viscosity very much 
larger than the actual, you can obtain instead of the discon- 
tinuity a front of a finite width, and as long as this width 
is small compared to the scale of the phenomenon, this 
should otherwise give no trouble??? Appendix A provides 
the original letter from Peierls. 


In 1948, Richtmyer^ said, “For strong shocks it has 
been customary to interrupt the normal calculating rou- 
tine at the discontinuity and perform a special calculation 
("shock-fitting") based on the Rankine-Hugoniot theory.” 
In that report, Richtmyer advocated for a new approach 
to simulate shocks using viscosity (real or fictitious) to 
smear the shock and achieve robust solutions without ex- 
plicitly fitting the shock. He said, "It was pointed out 
by von Neumann that what is clearly needed is to take 
into account in the equations the dissipative mechanisms 
that operate in a physical shock to convert mechanical en- 
ergy irreversibly into heat. This is the basis of the second 
proposal." With this approach, "shocks are automatically 
taken care of by ... introducing a (real or fictitious) dissipa- 
tion term...” Richtmyer” provided an analysis to show the 
Rankine-Hugoniot jump conditions were satisfied using 
an artificial viscosity. Using artificial viscosity had the 
advantage of simplicity and could be readily implemented 
on the computers of that time. 


By 1948, Richtmyer had developed a shock captur- 
ing scheme that incorporated artificial viscosity.^^?^?? 
Richtmyer added an artificial viscosity term to Peierls and 
Skryme's Lagrangian hydrodynamic scheme (without the 
shock-fitting boundary conditions) that was in terms of 


the initial coordinates. He also proposed a new scheme 
to solve the internal energy evolution equation, which in- 
cluded evaluating the EOS in terms of specific volume 
and temperature rather than entropy,” p = p(v, T). He 
also presented an approach to solve the entropy evolution 
equation. The methods by Richtmyer to evolve internal 
energy and entropy, along with some of his derivations,” 
were not discussed in the paper that von Neumann and 
he later published in 1950;”* likewise, that paper did not 
discuss the earlier work on hydrodynamic methods by von 
Neumann, Peierls, and Skyrme. We present the details on 
Richtmyer's 1948 scheme in Section VI.. 

The von Neumann and Richtmyer artificial viscosity 
scheme published in 1950% used a new staggered time 
integration method with a staggered spatial discretization 
that was in terms of the initial coordinates, where the 
latter came from the hydrodynamic scheme by Peierls 
and Skyrme, and the spatially staggered-grid approach 
came from earlier work by von Neumann.!® The origins 
of many modern Lagrangian hydrodynamic schemes can 
be traced to the pioneering work by Richtmyer and von 
Neumann, and the impact of that work will be discussed 
further in the Conclusions section. 


ILA. Historical documentation 


Thoroughly understanding the specifics of the early La- 
grangian hydrodynamic methods is somewhat challenging 
for a range of reasons. There are a limited set of histor- 
ical reports that document the work, and these reports 
were not necessarily intended to educate an outside audi- 
ence. The authors were inventing the syntax to describe 
the numerical approaches, so at times, their syntax differs 
from what is used in many journal papers and textbooks; 
likewise, many descriptions in the report centered around 
using the IBM PCAM e.g., punch cards, which is quite for- 
eign to modern-day readers. The documents also contain 
typographical errors. The Los Alamos staff at the time 
were experimenting with various methods of handling the 
shocks so many options were reported. They were also 
pioneering the transition from classical analytic methods 
to numeric methods, so many of the reports discuss both. 
The rest of this paper seeks to decipher this complicated 
story of the first numerical methods to simulate shock 
dynamics. 


IIl. JOHN VON NEUMANN'S METHODS 


von Neumann derived the Lagrange equations in terms 
of mass coordinates for the case of compressible, non- 
viscous, non-conductive hydrodynamics in his March 
1944 paper. In that paper, the position at time t is 
x = x(a, t) where “the substance contained in the interval 
a, a + da has the mass da". For the rest of the paper, 
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we will use the syntax by Skyrme, Peierls, and Richtmyer 
for consistency purposes, which deviates from what was 
used by von Neumann in his paper.'® Using the alternate 
syntax, the Lagrange evolution equations for a position X 
in 1D Cartesian coordinates is 


OX  ðp 
E 7 (1) 


| üm 
where the pressure is p and the mass is m. The current 
position is X, which deviates from today's commonly 
used convention that denotes variables in the initial co- 
ordinates with capital letters, and the current coordinates 
with lower-case letters. The specific volume in 1D Carte- 
sian coordinates is 

pus OX 1 


=F = (2) 


von Neumann presented a finite difference approach 
(based on central differences) to solve these governing 
equations numerically on a staggered-grid. This is the 
origin of the staggered-grid approach to solve the La- 
grange hydrodynamics equations that is still in use today. 
The focus of von Neumann’s original paper!® was on 1D 
Cartesian coordinates, but at the end of his paper, a short 
discussion was given on extending the scheme to spherical 
coordinates. 

We will now shift the focus to von Neumann’s work 
to extend his original scheme to spherical coordinates 
that was documented in laboratory reports such as the 
one by Skyrme.'? Shortly after arriving at Los Alamos, 
Skyrme adapted his shock-fitting treatment, mentioned 
in Peierls’ March 15, 1944 letter,” to work with von 
Neumann's hydrodynamic method, ? and he documented 
the method in LAMS-562.!? In what follows, we will 
adopt the notation of LAMS-562.!° 

The Lagrange evolution equation for the radius R in ID 
spherical coordinates is, 


= agi P (3) 
m 


Here m is the mass coordinate of an individual Lagrangian 
mass-point, “equal to (37) times the mass enclosed within 
the spherical shell corresponding to any point",'? The 
position of a mass element at a particular time, t, is defined 
as R. The specific volume, v, is defined as 
ðR? 1 
Mm F (4) 
von Neumann presented a staggered-grid finite differ- 
ence method to solve the governing evolution equations. If 
the time cycle is labeled by a superscript n, and the mass- 
points are labeled by a subscript i, the discrete equation to 
be solved is 


Figure 1: The staggered mesh used with von Neumann’s first 
Lagrangian hydrodynamic method 
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See Figure 1 for a graphical illustration of the variable 
locations. 
The time difference is 


At - De tlm Dei (6) 
The mass difference is 
E Tp] — mi 
Am Mipi mia 5 (7) 


A key assumption that von Neumann made was that en- 
tropy was constant everywhere away from the shocks. The 
pressure is determined using the equation-of-state from 
the entropy S and the specific volume 


Pi = pv; i$) (8) 


The discrete specific volume is 


yn, — GRUT = OT? 
one mpi — mi 


(9) 


Because the pressure could be calculated from the entropy 
and specific volume there was no need to solve the energy 
equation. The pressure was calculated at the mid-points 
between the mass-points, this was a staggered-grid La- 
grangian formulation. 

von Neumann pointed out the Courant time stability 
condition.!®? Peierls had empirically found the same 
condition.!! Skyrme terms it the ’L-condition’ in his 
report.'? Defining the local sound speed, 


fm 
Op|s 


were the partial derivative is taken with entropy S remain- 
ing constant, and should not be confused with the subscript 
s that is used later in the paper to denote a quantity at a 
shock. The Courant limit is 


(10) 


L= 


E u Me 


AR Am (11) 


|«: 
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A value of L = 1/2 is reported to be used.? As noted 
above, the equation of state (EOS) is needed to solve the 
difference equation. The final form of the EOS was deter- 
mined so as to limit the number of punch cards required, 
from LA-1058,”° 


plv, S) = Fi? +b (12) 
a a 


Here F and F are constants for a material, v is the spe- 
cific volume, and the coefficients a and b are known func- 
tions of S, the entropy. The entropy was calculated from 
the specific volume.? One thousand punch cards”° were 
used to tabulate all possible values of v/a and b. This 
completes the hydrodynamics equations that were used 
with von Neumann’s 1D spherical coordinate method. 

von Neumann proposed using this method for treating 
shocks by just letting the hydro deal with it, i.e., “com- 
pletely ignoring the possibility of shocks".!? He visual- 
ized the problem as a set of mass points connected by 
springs. He believed that the above hydrodynamic equa- 
tions would allow good approximations of a shock, but 
with oscillations behind the shock. “As soon as a shock 
has been crossed such oscillations must develop, and they 
have a perfectly good physical significance. They repre- 
sent the thermic agitation caused by the degradation of 
energy through the shock".? Only the average hydro- 
dynamic values were valid behind the shock. The test 
problems von Neumann carried out on the BRL PCAM 
in early 1944 were air blasts that had weak shocks.!5 The 
approach of letting an inviscid hydrodynamic method cal- 
culate the shocks, especially weak shocks will be referred 
to as the von Neumann method in this paper. When tried 
on strong shocks, von Neumann's method created more 
violent oscillations than expected. Peierls showed that von 
Neumann's method was only valid for weak shocks. ?:?? 
The conclusion was that von Neumann's method could 
not be used for strong shocks, and lead to them using 
shock-fitting methods. 


IV. RUDOLF PEIERLS' SHOCK-FITTING 


METHOD 


Peierls proposed an approach in his March 28, 1944 
letter to von Neumann that treated the "shock wave classi- 
cally".!?:2? This method required extrapolating the shock 
position to the next time step.’ Metropolis worked out the 
numerics of Peierls shock method in April 1944. Skyrme 
described Peierls’ shock fitting method in LAMS-562.'? 
Skyrme stated that the “main disadvantage of this method, 
apart from the fact that it introduced undesirable fluctua- 
tions behind the shock, was that a cycle (n+1) could not 
be completed until the positions of the points (i-1) etc. 
behind the shock had been found by numerical integra- 
tion. This introduced a delay between cycles".'? This is 


the same weakness as the extrapolation method Peierls 
described in his March 28, 1944 letter.” 

For strong shock problems, the pressure found in the 
“classical calculation" with Peierls’ method agreed closely 
with the average pressure calculated using von Neumann's 
method. However, the classical method showed that the 
high pressures from von Neumann's method on strong 
shock problems were spurious. 


V. TONY SKYRME'S SHOCK-FITTING METHOD 


Neither von Neumann's nor Peierls’ shock methods 
were really satisfactory. Peierls' classical method re- 
quired information ahead of the shock, and von Neu- 
mann’s method was only valid for weak shocks.!?^?? This, 
coupled with the observation of the constant in space pres- 
sure behind the shock, lead to the use of a third method for 
handling strong shocks that also reduced the “undesirable 
fluctuations behind the shock". 

Skyrme dramatically increased the sophistication of 
handling shocks with the intent of doing the best possible 
simulation. The shock-fitting method is documented in 
LAMS-562'? and in LA-1058.6 LAMS-562 was pub- 
lished in 1947 after Skyrme left Los Alamos, but it ap- 
pears to have been written in December 1944 or January 
1945. Similarly, LA-1058 was published in June 1949, but 
it contains a chapter written by Skyrme that thoroughly 
described his shock-fitting hydrodynamic method. 

In both LAMS-562 and LA-1058, Skyrme described a 
novel iterative approach to impose the Rankine-Hugoniot 
jump conditions at the shock front within a Lagrangian 
hydrodynamic calculation. Skyrme also presented a new 
Lagrangian staggered-grid hydrodynamic method based 
on the initial coordinates that is also referred to in an 
early letter by Peierls in 1944.7? Skyrme's hydrodynamic 
scheme was applied everywhere on the mesh in LA-1058, 
whereas in the case of LAMS-562, his hydrodynamic 
method was only used near the shock and von Neumann's 
method was used everywhere else. '? 

For strong shocks, Peierls’ shock fitting method was 
replaced with the shock fitting approach of Skyrme.'? 
Skyrme's method shown in LAMS-562P? has 63 steps! 
The "procedure used in computing the position of the 
shock is so complicated that it is impractical to use the 
I.B.M. machines for it, but it can be done satisfactorily 
with an ordinary calculating machine"? Even so, the 
shock procedure on the Marchant for a time step could be 
completed in about the same time it took the IBM PCAM 
to calculate a single hydro cycle.?® 40 

By March 1945 a standard shock procedure was in 
place, Skyrme's shock-fitting approach, combined with 
von Neumann's hydrodynamic method, was used for the 
strong shocks. '? The weak shocks were treated using only 
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von Neumann's method.'? By the late 1940’s, researchers 
had switched to using Skyrme's Lagrangian hydrodynamic 
method (which is in terms of the initial coordinates) in 
place of von Neumann's method (which used mass coor- 
dinates) but still with the shock-fitting approach. 

The remainder of this section will focus on describing 
the key equations used with the shock-fitting scheme by 
Skryme. We will present the key parts for the scheme 
documented in LA-1058 as it is the most complete. There 
appears to be only one minor change between the scheme 
documented in LAMS-562 and LA-1058, which will be 
highlighted. The method in LA-1058 assumes an initial 
mesh that is at rest with a single outward propagating 
shock starting at the origin (akin to the Sedov blast wave 
problem*':4?; the presentation that follows adheres to this 
assumption. As a note of caution to the reader, the original 
manuscripts by Skyrme contained typographical errors 
such as placing a minus sign on an equation, omitting a 
superscript on a variable, or not including a square root 
sign on a term. As such, the equations presented in this 
paper will in a few cases deviate slightly from what was 
printed in the original manuscripts. 


V.A. Governing 
equations 


The shock fitting hydrodynamic methods by both 
Peierls and Skyrme solved the governing analytic La- 
grangian hydrodynamic equations in terms of the position 
in the initial spatial coordinates r, which are termed total 
Lagrangian methods in modern literature. The position 
in current spatial coordinates is R. We follow the naming 
convention used by Peierls and Skyrme. The specific vol- 
ume of the material in the initial coordinates is vo and the 
pressure is p. The evolution of a position in the current 
spatial coordinates is, 


Lagrangian hydrodynamic 


aR R? dp 


ws coe (13) 
The specific volume at time t is 
om 
v= vo BS (14) 


It is reasonable to assume that the flow is adiabatic 1.e., 
there is no time for heat to transfer away. A continuous 
flow will be isentropic. The only change in entropy is 
across the shock; as such, the rate of change of the entropy 
with respect to time away from a shock is zero. 


V.B. Rankine-Hugoniot jump relations 


This subsection presents the Rankine-Hugoniot rela- 
tions used in the scheme by Skyrme.!”-*° The velocity of 
the shock is 


drs Ps — po 
dt vo — Us 


(5) 


The pressure at the shock is ps and the specific volume of 
the shock is vs. The shock pressure is a function of the 
specific volume and entropy, ps = (vs, Ss). The material 
velocity at the shock is 


OR 
u= (2E) = Voo) 19 
The acceleration of the shock is 
rs w dus (w? w?) (17) 
dt? 2us dt 
where 
ums ad Way — wg 
Ov VQ — Us 


The acoustic impedance is pc = 4/ — Sp g and is a material 
dependent variable that was tabulated. The rate of change 
of the specific volume at the shock is 


dus 
dt 


2W 2usvsW , 2 (Ov Op 
me d (a) tela) 


where Rs is the spherical radius at the shock front in the 
current coordinates. This equation is a function of the 
state of the material that is being shocked (po, vo) and a 
function of vs. As a reminder, all presented equations 
assume a shock traveling in the positive radial direction 
into a medium that is at rest. 


V.C. Finite difference solver for shocks 


The specific volume at time n is calculated using, 
3 


n 3 
aa (Ria) 
ma pp 
This equation is the identical to the strong mass conser- 
vation used in modern methods. The nodal positions in 
the initial coordinates are r; and the nodal positions in 
the current coordinates are R;, see Figure 2 for a diagram 
illustrating the mesh nomenclature. The finite difference 
equation for evolving nodal positions in the current coor- 
dinates away from the shock is 


(20) 
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The pressure p? i is calculated using an equation of state 
that is a function of [A 1 and entropy S, which is constant 
away from a shock. See Figure 3 for a diagram showing 
the variable locations with a notional shock profile. 

The finite difference equation shown in Eq. 21 will 
give rise to large oscillations near strong shocks. To cor- 
rect these oscillations, Skyrme proposed a shock-fitting 
approach that introduces a discontinuity with boundary 
conditions to the finite difference equation. At a shock, the 
pressure gradient behind the shock is based on the shock 
pressure instead of the neighboring pressure at p; , 1. 


41 -1 2 vo( RP? [ P5 — Pi- 
R =2RR -R — (At) Wm 
E rear 

(22) 


As mentioned earlier, the shock is moving in the positive 
direction relative to node R;. The shock pressure ps is a 
function of v? through an equation of state. The shock 
pressure and the specific volume for the shock are un- 
knowns and must be calculated using an iterative method 
and the Rankine-Hugoniot jump equations. The entropy 
only changes across the shock, so the flow is isentropic 
away from the shock. 


V.D. Solving the jump equations 


This subsection discusses how to calculate v7, which 
in turn is used to calculate pẹ through an equation of state. 
The shock pressure p? is used to evolve the position of 
the nodes R; in the current coordinate system and the 
shock position in the initial coordinates rs. The goal is to 
calculate the specific volume of the shock v? iteratively 
using 


n. ,n-1 At dus " dvs n=l 
OEO 


where e Js see Eq. 19, is a function of the unshocked 


state (po, vd), a function of the current shock position r2 
in the initial coordinates, a function of v7, and a func- 


TL TL 
tion of spatial derivatives at the shock (3) : and (32) i 
These latter spatial derivatives are calculated using, 


Ov (Us ~ Vi- i Op\ — 
ineo EUM 
Eq. 23 is iteratively solved until vý on the left side is equal 


TL 
to the v? that is used to calculate (5e ) on the right side. 


The shock position at the next time step is calculated 
using, 


Ps — Dii 


Pg. 


(24) 


Nie t 


R(r, t) 
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Figure 2: Skyrme's shock-fitting hydrodynamic method solved the governing equations in the initial coordinate system. 


Figure3: A notional shock profile is shown with corresponding 
variables used with Skyrme's shock-fitting method 


n+1 n At 
Ts = Tet a 


(25) 
The first derivative of the shock position with respect 


n 
drs 


to time ( di ) , see Eq. 15, is a function of vs. The 


second derivative of the shock position with respect to 


n 
time (4) , see Eq. 17 is a function of vs and Hos 
Skyrme found that calculating r7! using 
E At? (di 
rt = are re + ( =) (26) 


generated oscillations, so he used Eq. 25. 

Skyrme’s shock-fitting approach described in LAMS- 
562 is very close to what was presented in this section. 
One minor difference worth noting between the two re- 
ports (LA-1058 and LA-562) is that the shock position is 
calculated in LAMS-562 neglecting the second derivative 
term in Eq. 25. 


VI. ROBERT RICHTMYER’S HYDRODYNAMIC 
SCHEME 


An interesting historical point is at the end of Peierls’ 
March 28, 1944 letter to von Neumann where he suggested 
that an artificial viscosity could be used to slightly smear 
the shock discontinuity, making it numerically tractable.” 
Richtmyer implemented this idea in 1948.73 He and 


n drs\" | (dr V" ^| AP (drs\" 
2 dt dt 2 Vd? 


von Neumann made this the standard method for handling 
shocks in Lagrangian hydrodynamic codes with their 1950 
journal paper.?^ This is known as a smeared shock method 
or shock capturing method. Use of artificial viscosity 
allowed the intrinsic treatment of strong shocks without 
undue smearing of the weak shocks. The details on Richt- 
myer's scheme will now be discussed. 

The staggered-grid hydrodynamic scheme by Richt- 
myer?” solves the governing equations in 1D Cartesian 
coordinates in terms of the initial spatial coordinates zx. 
His approach builds on the one by Skyrme (subsection 
V.A.), and includes an artificial viscosity term and solves 
the specific internal energy evolution equation. This is 
the first recorded instance of solving the specific internal 
energy evolution equation with a hydrodynamics method, 
which is now the standard practice. The pressure is calcu- 
lated using the specific volume v and temperature T', so 
p = p(v, T). The prior hydrodynamic methods calculated 
the pressure in terms of the entropy S. The evolution equa- 
tions for the position in the current coordinate system X 
(in 1D Cartesian coordinates) and specific internal energy 
c are as follows, 


1dx Op 8v | 8v 

vo dt? Ox Mas Ot | OxOt (en 
de ðv | |Ov 3 (28) 
dt Pot at 


The artificial viscosity coefficient a is derived based on 
studying the Rankine-Hugoniot jump relations. 


1d 
PER Au 2 
2 UGUs 


(29) 


where ls is the shock thickness and vs is the specific vol- 
ume of shocked state. As an approximation, Richtmyer 
proposed to use the non-shocked value of specific volume 
vo in place of vs. The specific volume v is calculated using 


Ox 
ce 


(30) 


These governing equations are discretized using a stag- 
gered arrangement and central difference approximations 
of the derivatives. 
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The pressure is calculated using p", ; = p(v?, 1, T7. 1). 

itg itg ict 
Richtmyer also presented an alternate entropy-based ap- 
proach that evolves specific entropy s in place of the spe- 
cific internal energy. 


n+l  ,n n+l  ,mn 


an jeg ee Uni 
e$ ^ AT At 

In 1953, a code was implemented on the IBM Model 
II CPC.? It employed both the method by von Neumann 
and Richtmyer that used an artificial viscosity for shocks, 
and also Skyrme's shock-fitting method, but some hand 
fitting was still required.^ By the mid-1950's, the shock 
fiting method of Skyrme was fully automated, and the 
successor code had the option to use either shock-fitting 
or smeared shock methods for simulating strong shocks. 


v 


(34) 


VII. CONCLUSION 


The intent of this paper was to inform the broader com- 
munity on the pioneering work by von Neumann, Peierls, 
Skyrme, and Richtmyer to create robust hydrodynamic 
methods to simulate shock drive flows. These researchers 
developed multiple Lagrangian hydrodynamic schemes 
between 1943 and 1948, and most of these schemes were 
never published beyond laboratory reports. These early 
foundational advances led to subsequent, important com- 
putational hydrodynamic methods created by researchers 
at Lawerence Livermore, Los Alamos, and Sandia Na- 
tional Laboratories, and at laboratories and universities in 
the UK, France, and other countries. 

We discussed the original, 1D Cartesian coordinate 
staggered-grid Lagrangian scheme by von Neumann that 


Submitted to ANS/NT (2021). LA-UR-21-20144. 


was published in 1944, but garnered little attention by 
the broader community after it was published. That 1D 
scheme used central difference approximations for the 
derivatives. Soon thereafter, von Neumann proposed a 1D 
spherical coordinate Lagrangian hydrodynamic method, 
based on his earlier 1D scheme, to simulate a radially- 
outward propagating shock wave (akin to the Sedov blast 
wave problem) on the BRL IBM PCAM machine. We be- 
lieve this was the first shock hydrodynamics code, but that 
scheme was prone to spurious oscillations near a strong 
shock, often with amplitudes twice the mean pressure. von 
Neumann's original method was only suitable for weak 
shocks. 


At Los Alamos, von Neumann's hydrodynamic method 
was implemented in a computer program by Frankel, Nel- 
son, and Livesey, who also directed a cadre of IBM PCAM 
operators. These efforts resulted in the Los Alamos IBM 
code, which would become the first hydrodynamics code 
to accurately and robustly simulating strong shocks. 


The spurious oscillations in von Neumann's method 
motivated Peierls and Skyrme to investigate shock-fitting 
methods. We presented the details on a novel, essen- 
tially unknown shock-fitting Lagrangian hydrodynamic 
method created by Skyrme in 1944. It explicitly repre- 
sented a shock discontinuity on a discrete mesh and cou- 
pled together two separate higher-order solutions with the 
Rankine-Hugoniot jump conditions. With this approach, 
higher-order accuracy is achieved without spurious oscil- 
lations on problems with discontinuities. The scheme by 
Skyrme was an improvement over an earlier shock-fitting 
scheme by Peierls. The challenge with Skyrme's shock- 
fiting method was that it required an iterative solve to 
impose the shock boundary conditions, so it was not well 
suited for the IBM PCAM. Instead, the shock fitting was 
performed by hand using Marchant calculators with the 
results inserted into the IBM hydrodynamic calculation at 
the appropriate position and time. The scheme by Skyrme 
was fully automated on a computer by the mid 1950's. 


The shock-fitting approaches by Peierls and Skyrme 
were used with a staggered-grid Lagrangian hydrody- 
namic method that was in terms of the initial coordi- 
nates. That formulation quickly replaced Von Neumann's 
mass coordinate formulation. Since that time, a range of 
Lagrangian hydrodynamic schemes have been proposed 
that use the initial coordinates (or fixed reference coor- 
dinates), examples include the discontinuous Galerkin 
hydrodynamic method, ^^! the finite element hydrody- 
namic method,” and the residual distribution hydrody- 
namic method.?*?^ Skyrme also introduced a strong mass 
conservation equation that has been used in many subse- 
quent schemes.*6.52.55.56 


After WWII, researchers at Los Alamos were motivated 


to create a robust Lagrangian hydrodynamic scheme to 
simulate strong shocks and that could be fully run on the 
computers at that time. The alternative to the Skyrme 
shock-fitting scheme was a shock capturing scheme that 
used artificial viscosity to smear the shock over several 
cells. The use of artificial viscosity for shock capturing 
was first suggested in a letter by Peierls to von Neumann 
on March 28, 1944.” Richtmyer developed the quadratic 
formulation of artificial viscosity in 1948, calling it a 
“fictitious” addition to the hydrodynamics equations.?^?? 
Richtmyer added the artificial viscosity to Peierls’ and 
Skyrme's finite difference scheme that was in terms of 
the initial coordinates. Likewise, he introduced the now 
standard practice of solving the internal energy evolution 
equation, which included evaluating the EOS as a function 
of specific volume and temperature," whereas prior works 
used specific volume and entropy. 


The artificial viscosity approach was introduced to the 
broader hydrodynamics community, using Peierls’ arti- 
ficial viscosity nomenclature, in the well known 1950 
paper by von Neumann and Richtmyer."^ That paper pre- 
sented a new staggered time integration method with a 
spatially staggered discretization that was in terms of the 
initial coordinates. Peierls and Skyrme were the first to 
create a Lagrangian hydrodynamic scheme in terms of 
the initial coordinates following the staggered-grid ap- 
proach from von Neumann in 1944.!5 The staggered-grid 
discretization combined with artificial viscosity is the ba- 
sis of most Lagrangian hydrodynamic schemes.?9 9? In 
1955 at Los Alamos, Rolf Landshoff added a linear term 
to the artificial viscosity, 9? 64 and Marshall Rosenbluth 
suggested turning off the artificial viscosity during expan- 
sion.? These additions to artificial viscosity were first 
publicly discussed in the 1957 edition of Richtmyer and 
Morton.° Peter Lax and Burt Wendroff put the shock cap- 
turing method into conservative form in the 19605,6970 
In 1964, Mark Wilkins published the artificial viscosity 
in the form that is most widely used.’! In the same vol- 
ume, William Schulz described a tensor viscosity for use 
with a two-dimensional Lagrangian code." The stabil- 
ity analysis techniques and convergence theory of partial 
differential equations were established by von Neumann, 
Richtmyer, and Lax.?^7?7^ See Maattsson and Rider for 
a review of the development of artificial viscosity.” 


The creation of shock capturing schemes facilitated the 
development of two-dimensional Lagrangian shock hy- 
drodynamics codes, ?:?»71.72.7677 and eventually three- 
dimensional codes. ^9? ^ Modern staggered-grid La- 
grangian and artificial viscosity methods have evolved 
from the pioneering work of the 1940's.5- 75-94 The first 
cell-centered Lagrangian hydrodynamic method was pro- 
posed by Sergei Godunov.***° Research at LANL on 
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cell-centered Lagrangian hydrodynamic methods in the 
early 1980’s resulted in the CAVEAT code,5^55 which was 
followed by a myriad of work on cell-centered schemes in 
the 2000's.5?? The cell-centered work has led to notable 
improvements in the staggered-grid methods. 90-105 

This history shows the vital contributions made to the 
hydrodynamic methods theory by von Neumann, Peierls, 
Skyrme, and Richtmyer. It certainly justifies Hans Bethe's 
statement that the “collaboration of the British Mission 
absolutely was essential”! and Norris Bradbury's state- 
ment that “the British Mission supplied the major portion 
of experience in the field of theoretical hydrodynamics... 
It might also be pointed out that the United States was 
largely lacking in personnel experienced in this field of 
classical physics”.'°° This paper sought to correct histori- 
cal omissions and to communicate the complicated story 
on the origins of Lagrangian hydrodynamic methods for 
simulating shock dynamics. 
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A ORIGINAL LETTERS 


This appendix provides the original letters sent by 
Rudolf Peierls to J.R. Oppenheimer and J. von Neumann 
in March 1944. These letters were declassified and ap- 
proved for unlimited release, LA-UR-20-28408. Peierls 
made his first visit to Los Alamos on February 8, 1944, 
and subsequently moved to Los Alamos in June 1944 to 
lead the implosion theory group. The letters are the result 
of the discussions during that first visit to Los Alamos. 
These letters were found in the archives of the Manhat- 
tan Project held at the Los Alamos National Laboratory 
National Security Research Center. 

The first four letters are to J.R. Oppenheimer. 
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* March 9, 1944 discusses neutronics and hydrody- 
namics. 

March 15, 1944 mainly discusses hydrodynamics, 
giving the shock wave equations worked out by 
T.H.R. Skyrme. Skyrme joined Los Alamos on July 
31, 1944 and was a key member of the hydrodynam- 
ics modeling effort. 

March 21, 1944 discuses neutronics by Davison. 
March 24, 1944 discuses neutronics by Wilson. 


The last letter is to J. von Neumann on March 28, 1944 
discussing Peierls’ method for handling shock waves. The 
last paragraph of this letter is remarkable because Peierls 
suggests to von Neumann to use artificial viscosity to 
smear out shocks in hydrodynamics calculations. 
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9th March, 1944. 


TELEPHONE 
HANOVER 2-2460 
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Dr. J. enheimer, UNCLASSIFIED 
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Dear Oppenheimer, P.M. Lame Er 
I can now give you answers to some of the questions 
that I put to Wilson as a result of our recent discussions. 


(a) Perturbation Theory for Small Deviations from Spherical Shape. 


Wilson has looked into this and is writing up the 
theory. I understand from Mr. Skyrme, who has just arrived 
here from Birmingham, that Wilson considers this a very easy 
problem and, in faet, was in doubt whether there was not some 
more subtle point in our minds when we made the request, since 
the method he suggests to solve the problem is fairly obvious. 
I hope to get an account of this method from him shortly and 
shall, of course, then forward it to you at once. We can then 

see whether this takes care of the problem Bethe Was anxious to 
get solved. 


(b) Numerical Method for the Shock Wave Problem. 


l The present state of affairs is as follows: Host of 
the time lately has been spent on attempting to solve the equations 
for small deviations from the Neuman case to a high order of 
approximation, but it turned out that the equations for all but 
the first order become unmanageable. We believe, however, that 
a solution up to first order should give a sufficiently good 
basis for starting numerical integration and this had just started 
again when Mr. Skyrme left. It was, however, then too early 
to say what further difficulties might be encountered. In the 
meantime work has been done on finding an analytical form for 
the equation of state of air in the range concerned. This work 
was almost completed when Skyrme left put as he was working on 
this himself his departure may have caused some delay. I assume 
that, regardless of the fate of the numerical integration, this 
work on the equation of state should be pushed ahead as it would 
also be of use for Bethe's approximate method of solution. It 
is difficult to judge the time scale for this integration process 
as it may be fairly rapid if there are no snags. Our interpretation 
is that it should easily be completed in a time of the order of 
six months if everything runs smoothly, but it is quite possible 
that further difficulties may cause an appreciable delay. 
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Dr. J. R. Oppenheimer i -2- 9th March, 1944, 


Ee aite utet cae ai a a e ee ae! 


Mr. Skyrme has developed an alternative form of 
the equations in which the variables are r and r/Y where Y is 
the shock wave radius. This leads to fairly simple equations 
and to boundary conditions which areeasier than the usual ones. 
Skyrme will write down the equations which I will then send on 
to you. I don't believe that this method will be of particular 
interest in connection with your integration problem but it 
might give some improvement for the latter stages after part 
of the material has started moving outwards and when there is a 
definite shock wave. 


I hope you have in the meantime received the latest 
reports from Wilson to which I referred in a previous letter. 
These were transmitted to your Washington Office on February 9th. 
You may know that Chadwick has a full list of British Reports 
handed over, together with the dates when they were passed on. 

He may be able to give you information that would help to chase 
up any particular reports. 


Yours sincerely, 


E uL 


R. Peierls 
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HANOVER 2-2460 


15th March, 1944. 


Dr. J. R, Oppenheimer, 
P. 0. Box 1663, | «ance Y 


Santa Fe, New Mexico. 


Dear Oppenheimer, 


l have received an advance copy of the February Progress 
Report of the Birmingham Group. This, of course, will eventually 
reach you through official channels but as the final copy of the 
January Report does not even seem to have reached this country, 
these are evidently very slow and I want to mention to you 2 items 
that may interest you. 


On the shock wave problem, numerical integration has been 
resumed and led at first to disappointing results as the accuracy 
obtained seemed poor. ‘They have stopped to investigate the reason 
for this and will report to us any further results. In the mean- 
time, work on the equation of state for air is continuing. . They 
have worked out the general theory of a hemisphere conpletely 
surrounded by an infinite scattering container and the results are 
now being evaluated numerically. I would welcome some indication: 
on whether you think this prohem at allof interest and whether it 
should be pressed with any urgency. 


Dirac’ reports that his calculations on efficiency have 
now been generalized ‘to inelude the scattering in the container 
and extensive numerical work is now starting on this. This report 
is dated about two weeks after I cabled to Dirac,warning him about 
the error which Bethe has found in his method and I am not quite 
clear whether in the numerical work now beginning this has been 
corrected for. In any case I think it would be very helpful indeed 
if we could before long have the notes which Bethe promised to 
write up on this point. 


I mentioned in my last letter the new form of the shock 
wave equations which was suggested by Skyrme. We have looked into 
the question of their suitability for numerical.work and come to 
the conclusion that for purposes of numerical equations they are 
completely equivalent to the oid ones. . It is, however, of some 
interest that, in the variablesSkyrme uses, the equations are rather 
simpler than 1 would have expected and I enclose a note on the 
form of the equations for what it is worth. 


Yours sincerely, 
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In the usual notation 
R is position of particle initially at £ at time t. 
P 8 are pressure and density 
8, is initial density outside disturbance 


The equations are then 


$RO K >N l (1) 
SU 72 dY 


R? aR (2) 


dux 
where n. ^ l y- p (3) 


Boundary conditions; when veo, R-R=0 (4) 
when R- r = Y(t) (5) 
-Yw (6) 

Y. (7) 

-ae/fg = Ve — | (8) 

arh Ue ES (9) 


In terms of Z & 


2. | ez ». y 2 (10) 
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R= Y F(z) \ 
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Our equations are now 
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H= F T (20) 
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Boundary conditions 
At origin R=R=0O 
or R20, ÐR o when y-o (23) 
Ox 
At shock-wave 
when x 20,221, F-1 (24) 
(6) becomes y7 P= y. YH 
or (25) 
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Dear Oppenheimer, 


I enclose a copy of some notes by Davison on the 
boundary condition in the case of 2 media of different mean 
free path in case they are of interest io your group, i 
ought to explain that this is unfinished work - Davison 
abandoned this work temporarily at the time I left England 
and I took with me some manuseript notes of his which have 
been typed out here, You will see that he does not arrive 
at any very definite answer, but nevertheless the methods 


he uses may be of interest, 
Yours sincerely, 
fc. b uL. 


Dr. R. Peierls 
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24th March, 1944. 


Dr. JR. Oppenheimer, 
P.O. Box 1539, 
Santa Fe, New "Mexico. 


Dear Oppenheiner, 
Tepe 
r 

I enclose a Tupiy by Wilson on the question of Beatty 
spherical bodies which he wrote up at my request, following my 
conversation with Bethe on this point. I see that he has used 
diffusion theory. I believe it is not surprising that with the 
"elementary boundary condition" he obtains a very poor result 
because the effect will largely depend on the undisturbed density 
near the surface, which with the elementary boundary condition 
1S zeros 


I had hoped Wilson's caleulations would include a 
general formula based on perturbation theory applied directly 
to the integral equation, which would be rather similar to Dirac's 
perturbation method. It is a fairly easy matter to do so, but 
I am reluetant to go into the details of this as this would involve 
first looking into the point that Bethe made, namely, that Dirac's 

g method is incomplete and has to be supplemented by further terms, 

This no doubt would apply also to the present case, 


ae ” The enclosed copy of Wilson's report is our only copy 
pie and I would be grateful if you could have it returned to me when 
in due course you receive an official copy from Washington. 


"o I find that we did not at first obey the letter of 
the regulations by which receipts were to be sent with all SECRET 
—— in letters. I would therefore be grateful if some time you could 
M check that you have received my letters dated 1 February, 
«un 2lst February, 9th Ma , Which were sent without receipts. 
*M ia " 4t 
(27 Yours sincerely, 
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rode MDRL. Quis Y — 
— —bDr-4. von reunann 
Institute of Advanced Studies 
Princeton, N.J. 
Copied from Los Alamos 2 A 
Nationoi Laboratory Archives 


my 
S 3 year von hcuann, 
oc 
E» Following our recent conversation I an writing to give 
5 N you the procedure used in our case to solve the boundary conditions 
> NS in the case of a shock wave advancing into a stationary gas. i 
e believe the uethod can Le applied without ;uch trouble to the ore 
8 3 ! eral case, 
à & 
I shall assume, for siiiplicity, that r, t are the in- 
ic endent variables, and that the intervals chosen are equidistant 
o0th in space and tine, ‘urtier I all assunse that, at a given 
tiue, the position o? the shock wave coincides with the end of one 
of the intervals which we ise, Then, at tiat instant, ve can 
describe our state of anowledze by diagram 1. The variables known 
it each point from previous integration are indicated. R is tne 
Penal position, p tne pressure, « the velocity. Note tuat « at 
;,U-€-ne! jd > stly because of ig 


value for tue position of t 
Then, since the medium to 


y 


e shock wave at 
1¢ right is 


4 
+ } 
w 5 


sne 2, 3 
{e+ c : i " t ni d - 
l'sti;roed, we "ave R s yn, Yor C a N > 2, T z Yasa 

lni o $c 1e i + $ s$ ^i mien aed os 4 + 7 " 

,.u1S 1S related to the shock wave velocity Vasi by 

t 


ev = Ys Ya (1) 


^h I ; 
correct to second order inclusive, e being the time interval. Next 
we can find the density at a point half-way between rm., and | "S 


which is denoted by r* in diagram 2; 
* zx oJ. 
Y FEY. eet n) 


The density is given by 


(2) 
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where f^ is atuospieric density.  Eguition (3) is again correct 
up to terus of second order in the nuuerator. 


Je eonpare this with tne density behind the shock-wave 
at t = n+ l which isf(n +1,9.,,). The difference between this 
and f * is a snall correction of first order and need therefore be 
known only to first order of accuracy. Now note tnat the line from 
(n + 2, n - 3) to (n, m-l) is parallel to, and twice as long as, 
the line from (n + 2, r*) to (n + ere j hence to first order w 


can write 
* 
P ~ EESTI Ina) F Prr nen Pisis] (4) 


The accuracy of this can, if necessary, be iuproved by 
computing also f(n + 2, m - 5) - fin, 1073) ete. and extrapolating 
towards the shock wave, 


The density at the shocx wave f(n + 1, J..) is a function 
of Via, by the usual boundary condition. Since by weans of (1) (3) 
and (4) we can express tnis density as a function of Vw 4, we obtain 


an equation which can be solved for . 
Vives DOE ARCHIVES 
Having obtained Vw, , and hence Ypey we can also 
obtain u and p at tue shock wave. By interpolation we obtain 
R(m, n + 2) and u(a, n + 1) which are needed for continuing the 
integration process. Frown the value of p at the boundary we can 


will fcllow in future. By interpolation we obtain tue entropy for 


ue point m - 1. The interpolation requires tne knowledge of the 
actual position V,,, of the shock wave at t = n * 1. This can be 
obtained with sufficient accuracy fron values for previous tiues, 
tobetner with the knowledge ofN, | ‘ 

Pie process then repeats, except that the shock wave will 
now coincide with tne end of an interval and hence sone sligat 
uodification of the iiethod of extrapolation nas to be used. 


If the shock wave advances into a uedium which is itself 
in wotinn, we have to know the values of the variables u, p, f ahead 
of the shock front, in order to define the boundary conditions. It 
is probable that no trouble will be caused if one obtains these 
quantities by extrapolation, refining the process where necessary. 
For example, it is quite convenient to get the density ahead of the 
shock wave by substantially the same process as described above for 


tne other side. 


It is necessary for this approximation to have some 
Knowledge of the point V... since one has to extrapolate to this 
point, but probably this point can be predicted with sufficient 
accuracy beforehand. Otherwise a correction way have to be applied 
after the point has been found. : 


I am afraid these notes are very sketchy, but my wain 


object is to show that, if one can afford to take reasonably small 
intervals so that terms of higher order are indeed small, there is 
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nothing very terrible about having to introduce a discontinuity. 


Since this can be done for the correct thermodynamic 
properties of the medium, I think it has an advantage over the 
process that you described. One question on that process. Behind 
the shock front, this will involve fairly high-períod oscillations. 
In order to avoid appreciable errors, one wust therefore use very 
small tise intervals. Did you find this in practice an inconvenience, 
and how would the loss of tine caused by this compare with the gain 
owing to avoiding the discontinuity? DOE ARCHIVES 


Incidentally, have you thought at all about the following 
alternative way of avoiding discontinuity. In actual fact, the 
shock front has a finite width because of the viscosity and thermal 
conductivity of the medium. But artificially, assuming a viscosity 
very wuch larger than the actual, you can obtain instead of the 
discontinuity a front of a finite width, and as long as this width 
is small compared to the scale of the phenomenon, this should otherwise 
give no trouble, You pay for this by increasing the order of the 
equation as regards space derivatives, but the higher derivatives 
will be negligible except near the shock front. This is a macroscopic 
way of avoiding the discontinuity as opposed to your atomistic one. 
It is more couplicated than yours but you can put in arbitrary heats 
and equation of state for the wediun, 


Yours sincerely, 


ce: Dr. Oppenheimer 
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